function solve = solve(j,fp,param)
alpha = param(1:11);
TFP = param(12:13) ;    
    
fp_parfor = struct();
fp_parfor.mean_skills    = fp.schools_distrib(j,1);  % # mean skills; 
fp_parfor.sd_skills      = fp.schools_distrib(j,2);  % # standard deviation skills; 
fp_parfor.n_simulation   = fp.schools_distrib(j,3);  % # children in the economy;

if fp.do_counter_type == 4 %No Inequality;
fp_parfor.mean_skills = fp.national_mean;
fp_parfor.sd_skills   = 0.1;
end

if fp.do_counter_type == 5 %No Between-Neighborhood Inequality;
fp_parfor.mean_skills = fp.national_mean;
fp_parfor.sd_skills   = fp.national_sd;
end

if fp.do_counter_type == 6 %No Within-Neighb.Inequality (mean skills are neghborhood-specific and unchanged);
fp_parfor.sd_skills   = 0.15;
end


h = 1; %solve for the equilibrium first;
fp_parfor.init_draws = fp_parfor.mean_skills + fp_parfor.sd_skills.*(fp.init_draws{j}(1:fp_parfor.n_simulation,1,h) - mean(fp.init_draws{j}(1:fp_parfor.n_simulation,1,h)))./std(fp.init_draws{j}(1:fp_parfor.n_simulation,1,h));    

if fp.do_counter_type == 7
% Truncate the distribution at lower 10th percentile
skills_children_below_10  = fp_parfor.init_draws( fp_parfor.init_draws<prctile( fp_parfor.init_draws, 10) );
skills_children_above_10  = fp_parfor.init_draws( fp_parfor.init_draws>=prctile( fp_parfor.init_draws, 10) );
fp_parfor.init_draws( fp_parfor.init_draws<prctile( fp_parfor.init_draws, 10) ) = datasample( skills_children_above_10 ,size(skills_children_below_10 ,1) ) ;
end

fp_parfor.peers_forward  = fp.peers_forward{j}(1:fp_parfor.n_simulation,1:fp_parfor.n_simulation,h,:);
fp_parfor.peers_backward = fp.peers_backward{j}(:,1:fp_parfor.n_simulation,:,:);    
fp_parfor.PS_forward     = fp.PS_forward{j}(1:fp_parfor.n_simulation,h,:);

epsilon = 5;
toll= 0.05;
epsilon_skills = 5;

record_epsilon=[3 epsilon];

i = 1;

while (epsilon > toll  && epsilon_skills>toll) %Fix point;
    
fp_parfor.i = i; %iteration;   
fp_parfor.j = j; %neighborhood;
fp_parfor.h = h; %simulation;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;    
%Simulate Forward;    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

if i == 1

    %first iteration: set the initial guess about the distribution of
    %skills;
    fp_parfor.distrib_PS = zeros( size( fp_parfor.init_draws , 1 ) ,fp.periods-1);     
    fp_parfor.inv0{1} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);
    fp_parfor.inv0{2} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);
    fp_parfor.inv0{3} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);
    
    fp_parfor.inv1{1} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);
    fp_parfor.inv1{2} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);
    fp_parfor.inv1{3} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);    
    
    distrib_skills1 = exp(fp_parfor.init_draws(:,1));            
    distrib_skills2 = technology(1,alpha, distrib_skills1 , fp.Max_Time, fp_parfor.distrib_PS(:,1) ,  distrib_skills1 , TFP) ;
    distrib_skills3 = technology(2,alpha, distrib_skills2 , fp.Max_Time, fp_parfor.distrib_PS(:,2) ,  distrib_skills2 , TFP)  ;
    distrib_skills4 = technology(3,alpha, distrib_skills3 , fp.Max_Time, fp_parfor.distrib_PS(:,3) ,  distrib_skills3 , TFP)  ;       
    
    fp_parfor.distrib_skills(:,1) = distrib_skills1;
    fp_parfor.distrib_skills(:,2) = distrib_skills2;
    fp_parfor.distrib_skills(:,3) = distrib_skills3;
    fp_parfor.distrib_skills(:,4) = distrib_skills4;
    

    fp_parfor.grid{1} = prctile( distrib_skills1 , fp.percentiles_kid ) ;
    fp_parfor.grid_peers{1} = fp_parfor.grid{1} ;    
    fp_parfor.grid{2} = prctile( distrib_skills2 , fp.percentiles_kid ) ;
    fp_parfor.grid_peers{2} = fp_parfor.grid{2} ;              
    fp_parfor.grid{3} =  prctile(distrib_skills3 , fp.percentiles_kid );
    fp_parfor.grid_peers{3} = fp_parfor.grid{3};

else

    %solve forward to simulate the distribution of skills and investments;
    sf = solve_forward(param,fp,fp_parfor,h);
        
    previous_iter_distrib_skills = fp_parfor.distrib_skills;    
    
    fp_parfor.distrib_skills = sf.distrib_skills;
    fp_parfor.distrib_PS = sf.distrib_PS;
    
    fp_parfor.grid{1} = linspace( prctile(exp(fp_parfor.init_draws(:,1)) , 5) , prctile(exp(fp_parfor.init_draws(:,1)), 95), length(fp.percentiles_kid) );
    fp_parfor.grid_peers{1} = fp_parfor.grid{1} ;
    fp_parfor.grid{2} = linspace( prctile( fp_parfor.distrib_skills(:,2) , 5) , prctile(fp_parfor.distrib_skills(:,2), 95 ), length(fp.percentiles_kid) );
    fp_parfor.grid_peers{2} = fp_parfor.grid{2} ;              
    fp_parfor.grid{3} = linspace( prctile( fp_parfor.distrib_skills(:,3), 5 ) , prctile(fp_parfor.distrib_skills(:,3), 95 ), length(fp.percentiles_kid) );
    fp_parfor.grid_peers{3} = fp_parfor.grid{3} ;    
 
   fp_parfor.inv0{1} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);
   fp_parfor.inv0{2} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);
   fp_parfor.inv0{3} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);
    
   fp_parfor.inv1{1} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);
   fp_parfor.inv1{2} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);
   fp_parfor.inv1{3} = linspace(fp.Min_Time,fp.Max_Time,fp.inv_points);      

    
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Solve Backward;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

sb= solve_backward(param, fp,fp_parfor) ;      

fp_parfor.Hc_tp1_guess_prime = sb.Hc_tp1_guess_prime;

if i > 1 
eps1=zeros(fp.periods-1,1);
eps2=zeros(fp.periods-2,1);
eps_skills=zeros(fp.periods-2,1);


for t = 1:fp.periods-1
            
  eps1(t) = max(abs(fp_parfor.estim_policy0{t}(:) - sb.estim_policy0{t}(:)));
  
  if t==3

  end
  
    if t < fp.periods-1
        
        eps2(t) = max(abs(fp_parfor.estim_policy1{t}(:) - sb.estim_policy1{t}(:)));
 
    end
    
   if t > 1  
   deciles_skills_previous_iter = quantile(log(previous_iter_distrib_skills(:,t)),[0.05:0.1:0.95]);
   deciles_skills_current_iter = quantile(log(sf.distrib_skills(:,t)),[0.05:0.1:0.95]);          
   
   eps_skills(t) =  max(abs(deciles_skills_previous_iter-deciles_skills_current_iter));
          
   end    
    
    
end

epsilon = max([eps1;eps2]);
epsilon_skills = max(eps_skills);
%keep track of epsilon records;
ind_rec = mod(i-1,2) + 1;
record_epsilon(ind_rec) = epsilon;


end


fp_parfor.estim_policy0 = sb.estim_policy0;
fp_parfor.estim_policy1 = sb.estim_policy1;

fp_parfor.estim_value0 = sb.estim_value0;
fp_parfor.estim_value1 = sb.estim_value1;

if i ==1
fprintf('Solving Neighborhood %d:', j)
end
fprintf('%d ', i)

i = i + 1;
%record_epsilon

end

fprintf('Solved Neighborhood %d \n ', j)

Value0(:,:,:,:,h) = sb.Value0_actual;
Value1(:,:,:,:,h) = sb.Value1_actual;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%Simulate once the fixed point is achieved and create the simulated 
%dataset;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

fprintf('Simulating Neighborhood %d:\n ', j )

for h = 1 : fp.n_rip_sim

if h ==fp.n_rip_sim
fprintf('Simulated Neighborhood %d \n ', j)
end

    
fp_parfor.init_draws     = fp_parfor.mean_skills + fp_parfor.sd_skills.*(fp.init_draws{j}(1:fp_parfor.n_simulation,1,h) - mean(fp.init_draws{j}(1:fp_parfor.n_simulation,1,h)))./std(fp.init_draws{j}(1:fp_parfor.n_simulation,1,h));

if fp.do_counter_type == 7
% Truncate the distribution at lower 10th percentile
skills_children_below_10  = fp_parfor.init_draws( fp_parfor.init_draws<prctile( fp_parfor.init_draws, 10) );
skills_children_above_10  = fp_parfor.init_draws( fp_parfor.init_draws>=prctile( fp_parfor.init_draws, 10) );
fp_parfor.init_draws( fp_parfor.init_draws<prctile( fp_parfor.init_draws, 10) ) = datasample( skills_children_above_10 ,size(skills_children_below_10 ,1) ) ;
end


fp_parfor.peers_forward  = fp.peers_forward{j}(1:fp_parfor.n_simulation,1:fp_parfor.n_simulation,h,:);
fp_parfor.PS_forward     = fp.PS_forward{j}(1:fp_parfor.n_simulation,h,:);    
    
fp_parfor.i = i;    
fp_parfor.j = j;
fp_parfor.h = h;

sf = solve_forward(param,fp,fp_parfor,h);

sim_data(:,:,h) = sf.sim_data;
sim_data_network(:,:,h) = sf.sim_data_network;

end

solve.Value0 = Value0;
solve.Value1 = Value1;
solve.sim_data = sim_data;
solve.sim_data_network = sim_data_network;